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Abstract. We study a non-ergodic transition in a many-body Langevin system. 
We first derive an equation for the two-point time correlation function of density 
fluctuations, ignoring the contributions of the third- and fourth-order cumulants. For 
this equation, with the average density fixed, we find that there is a critical temperature 
at which the qualitative nature of the trajectories around the trivial solution changes. 
Using a method of dynamical system reduction around the critical temperature, we 
simplify the equation for the time correlation function into a two-dimensional ordinary 
differential equation. Analyzing this differential equation, we demonstrate that a 
non-ergodic transition occurs at some temperature slightly higher than the critical 
temperature. 
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1. Introduction 

Glassy systems exhibit distinctive phenomena, such as divergence of the viscosity pQ, 
a history-dependent response [2|, and an extreme slowing down of the time correlation 
of density fluctuations These phenomena have been studied extensively through 

laboratory experiments, numerical experiments, and theoretical analyses of models [Sj. 
Although each of these phenomena has gradually come to be understood individually, 
we are still lacking a unified picture that can account for all of them. In contrast to the 
rich variety of physical phenomena displayed by these systems, the types of theoretical 
methods that have been used in their investigation are few. Particularly, many analyses 
are based on the mode coupling theory (MCT) However, the MCT is not justified 
completely from a microscopic viewpoint, and also it is not regarded as the first-order 
description of a systematic formulation. Therefore, it is important to formulate a new 
method for analyzing the phenomena of glassy systems. 

Recalling the historical development of the theory of glassy systems, we realize that 
the MCT has been the most commonly used theory jB] because it describes various 
phenomena on the basis of the singular behavior of the two-point time correlation 
function of density fluctuations. This time correlation function exhibits a two-step 
relaxation as a function of the time difference and a plateau regime appears between 
the two steps in the relaxation. The MCT predicts the transition temperature for a fixed 
density (or the transition density for a fixed temperature) at which the plateau regime 
extends to infinity. This phenomenon is called a non-ergodic transition. Although 
the singularity might be an artifact of the approximation employed in the theory 
jO], there are evidences suggesting that the non-ergodic transition has been observed 
experimentally and numerically with fairly good quantitative correspondence to the 
theoretical prediction 0j. Noting such evidence, in the present paper, we theoretically 
study the behavior of the time correlation of density fluctuations. 

Let us briefly review the essence of the MCT. This theory provides a self-consistent 
integral equation of the time correlation function and the response function, ignoring 
vertex corrections jSj- Then, by analyzing this equation numerically, one can find a 
non-ergodic transition. Here, it should be noted that the memory contribution to the 
equation for the time correlation function is expressed in terms of the time correlation 
function. It is believed that such non-linear memory plays a key role in the non-ergodic 
transition. 

Now, it is interesting that behavior similar to that exhibited in the non-ergodic 
transition seen in glassy systems is easily obtained in ordinary differential equations 
possessing no nonlinear memory. As an example, let us consider Newton's equation of 
motion for a point particle moving in a one-dimensional space. Let x(t) be its position 
and V(x) be a potential function given by V(x) = —x 2 [(x — x*) 2 + e], with e > 0. We 
consider the particle motion under the condition that x — > as t — > oo. With this 
condition, it is easily found that when e is small, the particle starting from x(0) (> a:*) 
climbs the potential, slowly passes the non-zero maximum position, and then converges 
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to x = 0. The graph of x as a function of t exhibits a two-step relaxation behavior with 
a long plateau. Furthermore, if e — 0, the particle converges to the position x* and 
cannot approach the origin. This qualitative change in the behavior of the system as a 
function of e occurring at e = is analogous to that exhibited by glassy systems at the 
non-ergodic transition. 

Clearly, this simple Newtonian system has no direct relation with glassy systems. 
Nevertheless, it might be possible that there exists a variable representing the time 
correlation function for a glassy system that obeys an equation similar to that in the 
example discussed above. Such a variable would provide a characterization of the non- 
ergodic transition in that glassy system. Beginning from this speculation, in this paper, 
we determine such a variable and derive its evolution equation from a many-body 
Langevin model, employing several assumptions. We demonstrate that the obtained 
equation indeed does exhibit a non-ergodic transition. 



2. Model 



We investigate a system consisting of N identical Brownian particles suspended in a 
liquid. We denote the volume of the system by V and the temperature by T, and 
consider the thermodynamic limit, represented by V — > oo and iV — > oo, with the density 
p = N/V fixed. Let r*j(t) be the position of the i-th particle, with i = 1, 2, • • - N. We 
assume that the i-th particle interacts with the j-th particle through some interaction 
potential C/(|r*j — The motion of each particle is described by a Langevin equation 
with a friction constant 7. Then, the evolution equation of the fine-grained density field, 



dp(r,t) 
dt 



- [ d 3 r'p(r',t)p(r,tWU(\r 
.7 J 



+ -Vp(r,t) + 

7 



'2T 

7 



p(r,t)S(r,t) 



where H(= (E. x ,E. y ,E, z )) is zero-mean Gaussian white noise that satisfies 
(S a (r, t)E p (r', s)) = 8 aP 5{t - s)5(r - r') . 



(1) 



(2) 



This evolution equation can be derived exactly from the Langevin equation describing 
the motion of iV particles [7] . Explicitly, we choose the form of the interaction potential 

as 

E 2kEa 



U(k) 



(3) 



= \k\. Throughout this paper, f(k) is used to represent the Fourier transform 
f(k) = J d 3 re lk ' r f(r). All quantities in this system are converted into 
dimensionless forms by setting Ai = E — 7 = 1. Further, for the specific example 
that we analyze numerically, we consider the case in which a = 8.0 and A 2 = Ai. With 



where k 
of f(r): 
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this choice, U(r) has a minimum value —1.18 at r = 0.09, while U(r) — > oo as r — > 
and U(r) — > as r — > oo. Explicitly one can obtain 

^)-£-^fi + ^ V rEi(r) : erEi( - r) . <4) 



47rr tt 2 \ 2 dr y 

where Ei(r) represents an integral exponential function defined as Ei(r) = J_ dte t /t. 
We also fix the density p as p = 4.436. The temperature T is regarded as a control 
parameter. 

For later convenience, we define the quantity ^(k, k';t, s) = 5p(k,t)6p(k' , s)/p 2 , 
where 5p(r,t) = p(r,t) — p. The statistical average of ^(k,k';t,s) provides the two- 
point correlation function C(k,t) according to the relation 

(tt(fc, k'; t, s)) = (2n) 3 C(k, t - s)S(k + k'), (5) 

where we have assumed that the statistical properties of density fluctuations are 
symmetric with respect to spatial translations, rotations, and temporal translations. 
Although the last assumption is not valid below the glass transition temperature, we 
focus on the higher temperature regime, in which stationarity holds. Now, following 
the motivation of this study, we wish to obtain a differential equation that determines 
the time dependence of the correlation function C(k,t). For this purpose, as a trial, let 
us first consider the quantity d 2 C(k,t)/dt 2 . This quantity can be expressed in terms 
of C(k,t) and a four-point correlation function. Then, writing this equality, we obtain 
an evolution equation for C(k,t) by making the simplest approximation of omitting 
the third- and fourth-order cumulant terms of this four-point correlation function. The 
equation obtained with this approximation represents a mean field theory in the sense 
that fluctuation effects of \I/ are ignored when we consider its average value. As far 
as we know, this naive approximation was not employed in previous studies. It might 
be important to find some connection with more traditional approximations such as 
a Kirkwood superposition approximation. With this approximation, we obtain the 
following evolution equation for C(k,t): 

^(M)_^ w+ nVc ( M) 



dt 2 \7 7 

n 2 f d 3 ki ~ 
+ tW J^C{k u t)C(\k- kl lt) 

U(h) 2 (k ■ k,) 2 + U(kx)U(\k - fc x |) (k ■ kx) (k-(k- fc x ))l , (6) 

with t > 0. Here, we impose the boundary conditions lim^oo C(k, t) = and 
d t C(k,t)\ t=0+ = —k 2 Tpf^, which can also be derived from (jTJ). 

Here, it is expected that C(k,t) decays rapidly to zero in the limit k — * oo 
and exhibits no singularity as a function of k for any t. Therefore, we can expand 
C(k,t) in terms of the set of Hermite functions <p n (k), with n — 0, 1, • • • as C(k,t) = 
Y^=o^ri{t)4'nik). Then, (0) can be rewritten as 

°° dC (t^) 

Grim -Ti = B n (t), (7) 

m=0 
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°° dB (t) °° °° 

2^ Gnm ^7 = C m (t)L nm + M nm iC m (t)C[(t). (8) 

m=0 m=0 m,Z=0 

Here, we have introduced the new variable B n (t) in order to obtain a first-order 
differential equation. The coefficients L nm and M nm i in © are determined by the 
system parameters. For example, L nm is calculated as 

L nm = ^ P dk<P n (k)<P m (k) (U{k) + ^ , (9) 

while G nm in (JJJ) and (jHJ) is given by 

f°° 1 

G nm = -J dk k2 + k2 (j> n {k)(j> m {k), (10) 

where we have introduced a cutoff parameter A in order to allow expansion of the 
inverse of the Laplacian in terms of the Hermite functions. We take the limit A — > 
after the transition temperature is calculated for a system with A fixed. In the actual 
computations, checking the temperature in the cases A 2 = 10~ 2 , 10~ 3 , and 10~ 4 , we 
found that A 2 = 10" 3 is a sufficiently small value. The numerical values appearing in 
the analysis hereafter were obtained using A 2 = 10~ 3 . 
Now we describe a further approximation. Setting 

u = (Co, Ci, • • • , Ck-1: Bq, ■ ■ ■ , B K _i), (11) 

we approximate (|7j) and (JHJ by the following ordinary differential equation: 

G^ = T,u + Q{u,u). (12) 

Here, the 2K x 2K matrices G and E are defined as G nm = Gx+n^+m = G nm and 
^n,K+m = $nm, ^K+n,m = L nm , for 1 < 71, m < K (The other components are zero). 
Also, is determined from M nm i and represents a map from two copies of u to a vector. 
We chose K = 100 for the analysis reported below. (We checked the K dependence of 
L nm , G nm and M nm i, and confirmed that K = 100 is sufficiently large.) 



3. Analysis 

We investigate the system behavior as we change the temperature T. We first focus 
on the behavior of solution trajectories of (fT2*|) near the origin, u = 0. When |tt| is 
sufficiently small, this behavior is approximated by the linear equation obtained by 
omitting the nonlinear term in p2jl . Then, calculating the eigenvalues and eigenvectors 
of E, we can classify the solution trajectories. In the high temperature limit, from 
the form of L nm , it is easily found that all trajectories either approach or move away 
from the origin exponentially, because the eigenvalues consist of K pairs of positive and 
negative numbers. We found in the numerical computations that when the temperature 
is decreasing there exists a temperature To below which one pair of eigenvalues becomes 
zero. (We have T =19.535 in the example we study.) Thus, the solution trajectories 
in the 2i^-dimensional phase space undergo qualitative change when the temperature 
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passes through T . The corresponding eigenvectors at T , which we denote $ o an d $01 ) 
play a prominent role in the subsequent analysis. 

It is known that this type of qualitative change near T = T can be described by 
a nonlinear differential equation for the amplitudes of these zero eigenvectors jH]- We 
now attempt to derive such an equation. Let £ be £ at T = T . Then, we find that 
$00 an d $01 satisfy the relations 

S $oo = 0, (13) 

E $oi = $oo- (14) 

Now, setting e = (T — T )/T , we consider the solution trajectories satisfying u — > in 
the limit t — > oo for the system with small positive e. For such trajectories, we expand 
u in eigenvectors of E as 

u{t) = Ai(t)e*oo + A 2 (t)e 3 / 2 $ i + r?(^i (*), (15) 

where r) represents the contribution to u that is not from the zero eigenvectors, and the 
amplitudes of the other eigenvectors quickly decay to the values determined by A\ and 
Ai. Then, the time evolution of the variables A\ and A2 is described by autonomous 
equations of the form 
d A (t) 

-^l = n t (A 1 (t),A 2 (t)). (16) 
Now, substituting (fTHjl into (fT2^1 and using expansions 

TJ = ^2^(3/2) + ^(2) + ^(5/2) + . . . ? (17) 

A, = eV 2 0f /2) + e 1 ^ + 6 3 / 2 fif /2) + • • • , (18) 

we can determine the functional forms of r}(Ai, A 2 ) and Qi(Ai, A 2 ) using a perturbative 
method. Note that the dependences on e in (fT5J) . (fTTjl . and (fT£J) are chosen so that the 
perturbative expansion can be carried out in a systematic manner. For example, all the 
terms proportional to e 3 / 2 in the expression obtained with the above procedure yield 

G^ 00 Q i ) l/2 \A 1 ,A 2 ) = A 2 ^ 00 + ± r 1 ^ 2 \A 1 ,A 2 ). (19) 

Because S has a zero eigenvalue, it is not invertible. Thus, when we regard (JT§j) 
as a linear equation for r]^ 3 ^ 2 \ there are either no solutions or an infinite number of 
solutions. In order to carry out the perturbative calculation systematically, we impose 
the solvability condition under which rj^ 3 ^ 2 ^ can be obtained uniquely. This solvability 
condition determines Q^ 2 ^. Applying similar procedures to each of the equations 
appearing successively, which are proportional to e 2 , e 5//2 , ■ ■ •, we derive , ^2^\ 
• in an iterative manner. Combining the contributions up to Vl^ and 
Vt^ 2 \ we obtain 

= aA 2 , (20) 

dt 

^ = Ml + b2 Al + 5 3y4 3 + ^ (21) 
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where each coefficient can be calculated as a function of T. For the example we analyzed 
numerically, we found a = -0.061186c 1 / 2 , b x = -1.108178e 3 / 2 , b 2 = -0.0005666 1 / 2 , 
b 3 = -0.011904e 3/2 and b 4 = 0.836102e 3/2 . (The details of the calculation will be 
reported in another paper 




A 1 A, 



Figure 1. Solution trajectories of (|2U[1 and (|21|) in (A\, A2) space for three values of 
e. (a) e = 2.520 x 1CT 3 , (b) e = 2.294 x 1CT 3 , and (c) e = 2.285 x 1(T 3 . The red curves 
represent physical solutions that satisfy lim^oo Ai = 0. Inset: A\ as a function of t. 
The red curves in the insets correspond to the red trajectories. 

Using this result, we plot trajectories of (Ai(t) , A 2 (t)) in Fig. 1 for three values 
of e. It is seen that for e = 2.520 x 10~ 3 , the trajectory approaches the origin with no 
particularly notable features. By contrast, for e = 2.294 x 10~ 3 , there appear two fixed 
points, and the trajectory approaches one of these. However, in this case, it eventually 
begins to move away from this point and converges to the origin in the t — > 00 limit. As 
shown in the inset of Fig. 1, the part of the trajectory near the fixed point corresponds to 
the plateau in the graph of Ai as a function of time. Finally, for e = 2.285 x 10 -4 , we find 
that the trajectory converges to the non-zero fixed point as t — > 00. The appearance of 
the connection between this new fixed point and the origin is called a saddle connection 
bifurcation in dynamical system theory. The inset of Fig. 1 shows that the saddle 
connetion bifurcation in the two-dimensional system studied here corresponds to the 
non-ergodic transition of the time correlation function. 

4. Concluding remark 

We have elucidated the nature of a non-ergodic transition using a method of dynamical 
system reduction to treat the time correlation function. The key idea employed in 
our analysis is to focus on the zero eigenvectors $00 an d <fr i at the temperature T . 
Using these eigenvectors, we can obtain a simplified description of the large- dimensional 
dynamical system (j!2)l in the form of the two-dimensional system represented by (|2T)|) 
and ()21|) near the temperature T . We find that the non-ergodic transition in the original 
system is manifested as the saddle connection bifurcation in the two-dimensional system. 
Before ending this paper, we make three important remarks. 

First, our theoretical result for the system with the potential (jSJ) has not been 
compared with numerical experiments yet. The comparison will be a future work. Here, 
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we note that our method can be applied to Langevin equations with other interaction 
potentials. As an example, for a model recently studied in Ref. JTU]; we confirmed the 
non-ergodic transition by both our method and a direct numerical experiment, though 
there is no precise quantitative correspondence between them. We also note that we do 
not understand the dependency of the results on potential forms. 

The second remark is that ((TJ takes the same form as the density evolution equation 
in the over-damped limit obtained with the dynamical density functional method, where 
the potential function in Q is replaced with the direct correlation function [Hj. It is 
thus expected that our method can also be applied to this type of density evolution, 
though the temperature dependence becomes complicated. 

Finally, recall that the equation we analyzed, ([12)1 . was obtained by ignoring 
fluctuations of ^(k, k'; t, s). Thus, it is important to extend our treatment to one that 
accounts for fluctuations, in which (j!2j) is derived as the zero-th order description. This 
problem is related to that of the growth of the length scale that appears near the glass 
transition point, which has been studied recently \1'2\. We are now in the process of 
developing such a theory. 
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